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ABSTRACT 

We report the discovery of a transiting, gas giant circumbinary planet orbiting the 
eclipsing binary KIC 4862625 and describe our independent discovery of the two tran- 
siting planets orbiting Kepler-47 (Orosz et al. 2012). We describe a simple and semi- 
automated procedure for identifying individual transits in light curves and present our 
follow-up measurements of the two circumbinary systems. For the KIC 4862625 sys- 
tem, the 0.49 ±0.018 rjupUer radius planet revolves every ~ 138 days and occults the 
1.14±O.14M0, 1.59±O.O6i?0 F8 IV subgiant primary star producing aperiodic tran- 
sits of variable durations commensurate with the configuration of the eclipsing binary 
star. Our best- fit model indicates the orbit has a semi-major axis of 0.56 AU and is 
slightly eccentric, e = 0. 1 . For the Kepler-47 system, we confirm the results of Orosz 
et al. (2012). Modulations in the radial velocity of KIC 4862625A are measured both 
spectroscopically and photometrically, i.e. via Doppler boosting, and produce similar 
results. 
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Introduction 



For decades the science fiction community has imagined that planets can orbit binary stars, 
yet only recently have such systems actually been detected. Timing variations either in the rotation 
period of a neutron star member of a binary system (Sigurdsson et al., 2003) or in the stellar 
occultations (when the two stars eclipse each other) of eight eclipsing binary (EB) systems (Deeg 
et al. 2008; Lee et al. 2009; Beuermann et al. 2010, 2011; Potter et al. 2011; Qian et al. 2011; 
Qian et al. 2012) have been interpreted as the gravitational perturbation of additional bodies on the 
binary stellar system, suggesting the presence of a total of 12 circumbinary (CB) planets on wide 
orbits with periods of tens of years. 

The lower limits on the masses of all twelve objects, however, fall in the super- Jupiter regime, 
making their planetary nature uncertain. Furthermore, the orbital stability of some of the multi- 
planet circumbinary systems (HW Vir, HU Aqr and NN Ser) have been studied recently, showing 



that so m e of them are on highly unstable orbits (IHomer et al.ll201 Ie iHinse et al.ll2012l ; iHomer et al 
2012alJbl : lOozdziewski et allboii beuermann et al.ll2012l) . 



Doyle et al. (2011) announced the first direct evidence of a Saturn-sized planet that tran- 
sits both members of an EB eclipsing binary, specifically Kepler- 16. Since then, five more CB 
planets have been announced: Saturn-sized planets that transit Kepler-34b and Kepler-35b (Welsh 
et al. 2012), the first Neptune-sized planet that transits Kepler-38 (Orosz et al. 2012b), and two 
Neptune-sized planets that transit Kepler-47 (Orosz et al. 2012). The latter system is the first binary 
discovered to have two planets, one Neptune-sized on a 300 day orbit and the other Earth-sized on 
a 49.5 day orbit. 

Substantial efforts in theoretical modeling indicate that planets such as these should not be 
uncommon. Simulations of dynamical stability show that beyond a critical distance, CB planets 
can have stable orbits in practically all binary configurations. The critical distance is on the order of 
a few binary separations (Dvorak 1986; Holman & Weigert 1999; SchoU et al. 2007; Haghighipour 
et al. 2010; Schwarz et al. 2011; Doolin & Blundell 2011). The orbits of Kepler- 16b, 34b, 35b 
and 38b are indeed outside the critical orbital semi-major axis, but only by 21%, 24%, 14% and 
26% respectively (Welsh et al. 2012; Orosz et al. 2012b). Kepler-47b, while notably farther from 
the instability region (Orosz et al. 2012), is still not too far out, suggesting such "reaching for the 
limit" behavior to be typical of CB planets. 



The fact that these four planets are so close to the theoretical limit for stability may suggest 
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that their host systems had an interesting dynamical history where migration and/or planet-planet 
scattering may have played a significant role in sculpting their present architecture. Formation 
and evolution theory of giant planets around binary stars has been studied extensively (Pierens and 
Nelson 2007, 2008a, 2008b), providing a number of outcomes that depend on initial conditions. 
Simulations of gas giants have shown that Saturn-size planets (like Kepler- 16b, 34b and 35b) 
stabilize at a 5: 1 orbital resonance and may be very common, compared to Jupiter- size planets that 
are either scattered out of the system or gradually drift outward into the disk. Single, Neptune-size 
planets (such as Kepler-38b) migrate and stop at a distance of about three times the binary stars 
separation, leading the authors to suggest that "the cavity edge of the precursor CB disk appears 
to be an excellent place to look for low mass planets in close binary systems." If there are two 
Neptune-size planets in the system, they become locked in a mean motion resonance, while a 
five-planet system is either disrupted or, in one simulation, also ends up in a resonance, implying 
that such multiple Neptune-size planets in resonant orbits may be indeed common. Models for 
the formation and evolution of terrestrial planets around binary stars (Quintana & Lissauer, 2006) 
have shown that, in the presence of Jupiter at 5 AU, CB terrestrial planets can readily form around 
a wide variety of binary systems. At least one terrestrial planet forms in all simulations presented 
by the authors. While the final masses of all simulated terrestrial planets vary little, the outcome for 
the architecture of the planetary system is very dependent on the parameters of the stellar binary. 
Highly eccentric binaries tend to harbor fewer, more diverse suite of planets compared to binaries 
with very low eccentricity, a prediction that can be addressed by the addition of more pictures to 
the family portrait of the five Kepler planets. 

More than 20 years ago Borucki and Summers (1984) proposed monitoring EB systems to 
search for planets because nearly edge-on inclination significantly increases the probability tran- 
sits. At the time is was not practical to monitor targets continuously over many days (Kepler- 16b, 
for example, has an orbital period of 230 days). On March 6 of 2009, 380 years after Johannes 
Kepler predicted the transit of Venus across the disk of the Sun, NASA launched the appropriately 
named Kepler Mission to search for Earth-like planets in the habitable zone of Solar-type stars and 
to determine their occurrence frequency (Borucki et al. 2010). To achieve this, a 0.95 Schmidt 
telescope on a Heliocentric Earth trailing orbit continuously and simultaneously monitors about 
150,000 stars in the visible range from 423 nm to 897 nm over a 100 square degrees patch of the 
sky in the Cygnus region where ~ 60% of the stars are G-type stars on or near the Main Sequence. 
Utilizing the transit method, the instrument searches for periodic dips in the brightness of a star 
caused by a planet transiting across its disk. The Kepler mission has been remarkably successful 
in finding transiting planets, discovering more than 2300 planet candidates (Borucki et al. 2011; 
Batalha et al. 2012), 77 of which have been confirmed by the time of writing. Amongst this trea- 
sure throve of data are also a set of 2165 EB systems (Slawson et al. 201 1), the main focus of our 
work. 
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Searching light curves of EB stars for transits of a third body is non-trivial. In addition to 
the significant limitations associated with intrinsic stellar variability and instrumental artifacts, CB 
planets have inconstant transit times, durations, and depths, all of which that depend on the phase 
of the binary system (Schneider & Chevreton 1990). To transit one of the stars, the planet must "hit 
a moving target" (Orosz et al. 2012). A benefit is that these transit signatures cannot be attributed 
to the stars themselves, to a background EB, or to other unrelated astrophysical or instrumental 
events, strongly supporting the CB-planet hypothesis. One challenge is that traditional transit 
searching algorithms, e.g. Box-fitting Least Squares (BLS) (Kovacs et al. 2002), used to detect 
periodic, box-like signals in the light curve of a single star are not optimized for finding transiting 
CB planets due to the unique nature of their signal. Several methods for the detection of transiting 
CB planets have been proposed. One approach is based on simulating light curves produced by 
an exhaustive search of possible orbits of CB-planets and fitting them to the data (Doyle et al. 
2000; Ofir 2008). Carter & Agol (2012) have developed the Quasi-periodic Automated Transit 
Search QATS algorithm, which is similar to BLS but optimized for aperiodic pulses. QATS has 
been successfully applied to CB-planets. Orosz et al. (2012) report that QATS failed to detect 
the outer planet Kepler-47c due to decreasing sensitivity for longer periods. While transiting gas 
giants cause a dimming of their host star large enough to be seen by eye in light curves, the transits 
of smaller planets can be easily missed by visual inspection. 

The initial discovery of CB transits together with the availability of exquisite Kepler data 
inspired us to develop a semi- automated procedure to identify aperiodic transits. We describe the 
procedure, which is based on the established BLS algorithm but modified and applied in a novel 
way. We applied it to finding transiting planets around EB stars listed in the Kepler catalog of 
Slawson et al. (2011). We examined the detached EB systems and identified several candidates 
that exhibited additional transit-like features in their light curves. Here we present the independent 
discoveries of two CB planets Kepler-47bc and KIC 4862625b. 

In this paper, we will describe the analysis as a linear, deductive process, although it is in- 
herently iterative, with one aspect feeding back into an earlier part. For brevity and clarity, we 
do not emphasize the iterations. This paper is organized as follows. In Section [Z2l we describe 
the procedure used to discover the two CB systems, followed by radial velocity measurements and 
spectra of the host stars in Section 13.21 The three-body dynamical model used to explore the pa- 
rameter space of the two systems is presented in Section |7]aka|71 We present our results in Section 
[6l discuss them in Section [8] and draw conclusions in Section [9l 
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2. Kepler Photometry 
2.1. Detrending 

We began with the long-cadence (~30 min) PDCSAP flux of KIC 4862625 generated by the 
Kepler mission for the publicly available^ quarters 1-9. Using the ephemeris of Prsa et al. (2010) 
and examining the phase-folded light curve, we flagged data points within 0.12 days of the centers 
of the primary or secondary eclipses, or within 0.5 days of the planetary transit events]^ There 
were 437, 387, and 238 points in each category, respectively, corresponding to 1.26%, 1.12%, 
0.69% of the total, 34681 data points. To remove the instrumental discontinuities in flux created 
by the quarterly rotation of the Kepler focal plane, we divide each quarter's data by its median. We 
flagged six points that differed from the rest by more than 0.2% in normalized flux. 

The unflagged flux is sinusoidal with a period of prot ~ 2.63 d, which we attribute to modu- 
lation caused by the rotation of star A. To detrend the flux variations attributed to rotation, for the 
purposes of providing eclipse and planetary transit light curves normalized to unity, we fit a unique 
sine wave in the local vicinity of each of the 34681 data points, in each case using unflagged data 
within Prot/2 of each data point. Each of the 34681 sine waves had a fixed period, prot, and we fit 
for three parameters: the mean, the amplitude, and the phase. We iterated the procedure, adjusting 
the value of prot, until the gradient in the phase shifts of the fitted sine waves over the entire data 
span was zero, indicating a best-fitted average rotational period, prot = 2.6382 ±0.0037 d, where 
the quoted "uncertainty" is the standard deviation of prot fit piece- wise over quarterly time spans. 
Finally, we divided each data point by its best-fitting sine wave evaluated at the time of each par- 
ticular data point. We use the resulting detrended, normalized light curve for subsequent analysis, 
and we use the mean values of the sine waves for the analysis of the Doppler boosting (Section 
14.41) . Exclusive of flagged data, the RMS of the residuals of the detrended, normalized light curve 
about unity is 222 ppm. 



2.2. Box-fitting Least Squares, BLS, for Single Transits 

We invented a semi-automated procedure to identify individual, aperiodic transit-like features 
in a light curve. The procedure automatically finds square-wave or "box-shaped" features within a 
light curve. 



http://archive.stsci.edu/kepler/ or http://exoplanetarchive.ipac. caltech.edu| . 



^This is one example of the iterative analysis: first we detrended the light curve, then we identified the planetary 
transits, then we detrended the light curve again with the transit points flagged. 
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Fig. 1. — Characteristics of the light curve of KIC 4862625. A 23-day portion of the Kepler light 
curve illustrates various phenomena. Instrumental effects are labeled above the data; astrophysical 
effects are labeled below the data. We flagged two data points as a "glitch," and a gap in the data 
is visible at BJD = 2455217. The eccentric orbit of the EB is apparent from the fact that primary 
eclipse is not centered in time between the two secondary eclipses. The primary eclipse is ~ 1.4% 
deep. The red line illustrates the sine-wave fit to the rotational modulation of the light curve. The 
secondary eclipses, the planetary transit, and the 2.6-day rotational modulation all have similar 
amplitudes, ~ 0.1%. 
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Using the raw EB light curves SAPFLUX from the Kepler database we first normalize them 
and remove the EE's eclipses using the Box-fitting Least Squares method (BLS, Kovacs et al. 
2002). Next we detrend the data, using only those points with a SAPQUALITY flag of 0. The 
detrending is non-trivial and has to be done on a target-by-target basis, as each binary star has 
complex baseline variability that spans timescales from hours to days. We use an iterative fit with 
a high-order Legendre polynomial on each quarter, further broken down into smaller segments 
according to the data flags. 

Within each segment of the light curve, the procedure automatically identifies the center, 
width, and depth of the two most-significant box-shaped features, one positive and one negative. 
The latter is a "transit" candidate and the former is an "anti-transit." Because transits are nega- 
tive features in the residuals of a detrended light curve, we can validate empirically the statistical 
significance of the transit candidates by comparison with the anti-transits extracted from the same 
data. Doing so helps us to not be overwhelmed with false-positives. 

The number M of light-curve segments is not particularly critical. The segment length should 
be longer than any actual transit. For CB planets, in principle the transits can be as long as half an 
orbit of the EB, in the case where the planet and the star are traveling in parallel at nearly the same 
projected rate (Schneider & Chevreton 1990). Such transits will be rare and generally accompanied 
by shorter transits, created when the planet and star are moving in opposite directions. The segment 
length should be shorter than the orbital period of the planet, lest only one of two genuine transits 
be identified in a segment. For CB planets of P-type (Dvorak 1982), the planet's orbital period 
must be longer than the period of the EB, and for orbital stability the CB planet's period must be at 
least a few times longer than the EB's period (Section|7]). From these two limits, the segment length 
should be between one half and a few times the orbital period of the EB. However, in practice, the 
intrinsic variability of the stars in the system prescribe the length of the segment. In any case, to 
prevent a transit being split by the segment boundaries, we analyze each light curve at least twice 
with the boundaries of the segments shifted. 

In all cases, detrending of the light curve is crucial to our method. Very short-period EB stars, 
contact, and semi-detached EBs are difficult to study with our method. They are highly variable 
and there are few measurements in between the stellar eclipses, forcing us to use segments much 
larger than the binary period. 

To avoid systematics effects that might mimic a transit, a merit criterion is necessary. As a 
convenience, one can use the BLS algorithm to find the individual transit- and anti-transit-features 
in segments of the detrended light curve's residuals. Because BLS is designed to find periodic box- 
like features, one option is to replicate N times the segment under study to form a periodic light 
curve. One can then use BLS to search for the most significant transit- and anti-transit features 
with that single period defined by the replication. The number N of replications is not particularly 
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important. To find the anti-transits, we simply invert the sign of the residuals of the light curve and 
search a second time. Although it would be more efficient to extract the most-significant positive 
and negative features in one pass, the execution time is trivial and enormously smaller than in the 
traditional BLS which must loop over thousands of trial periods, and smaller than the detrending 
also. Alternatively one can modify the original BLS computer code to analyze a single segment 
of unreplicated data; we have implemented each variant. The only potential difference is a tuning 
parameter within BLS, namely the minimum number of data points within an acceptable box-like 
feature, which must be adjusted commensurate with the number of replications N. 

After the automated procedure identifies the M most-significant pairs of transit- and anti- 
transit-candidates from the M segments, we examine the ensemble for outlier transit-candidates in 
a manner similar to that described by Burke et al. (2006). Burke et al. validated transit-like features 
with an ensemble of features reported by BLS for hundreds of stars observed simultaneously. 
Because of the large number of observations in each star's Kepler light curve, we validate transit- 
like features in each light curve with the ensemble of transit- and anti-transit features from only 
that particular light curve. 

For completeness, we briefly describe Burke's method here. Assuming dimming features 
(transits) and brightening features ("anti-transits") are due to systematic effects, it is reasonable to 
expect that there will not be a strong tendency for such effects to produce dips versus blips. In other 
words, typically ^{ji^transit '^i^^ ^e similar to ^(y})anti-transit s^ch segment. On the other hand, 
a highly significant transit signal is an outlier in a ^{jiyransit ^'^d ^{x^)antitransit diagram (Figure [21). 
The segments where noise dominates (black crosses) cluster in a cloud of points with similar values 
for ^(^yi^ransit ^{y})anti-transit^ ^ut the Segments Containing the tertiary transits (red crosses) are 
well separated. As seen from the figure (and depending on the merit criterion) there can be a 
significant number of outliers, requiring a human eye to check the segments which triggered the 
routine. The number of triggers are highly dependent on the detrending of the baseline - quiet stars 
(like Kepler- 16) have very few outlying points (planetary transits) where smaller transit signals in 
more variable binaries (Kepler-47) will be accompanied by a greater number of false positives. 

Figures |2l [3] and |4] show examples from the light curves for KIC 4862625, Kepler-47 and 
Kepler- 16. The three figures were presented on February 8, 2012 to a committee of Johns Hopkins 
University as part of the graduate student matriculation procedures for one of us (V.B.K.). As such, 
it used all the Kepler data that was publically-available at that time, namely Quarters 1 to 6. 

To evaluate the sensitivity of our method, we insert fake transits in the light curves and study 
their recovery rate as a function of the transit depth. An example for Kepler- 16 is shown on 
Figure |4] where the black crosses represents segments with zero fake transits, blue diamonds are 
segments dominated by systematic effects, red symbols are due to Kepler- 16b and the yellow stars 
correspond to segments in which fake transits with a depth of 240 ppm (Earth-like signal) were 
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superposed. The fake transits have variable depths, durations and period to simulate that expected 
for a CB planet. Also they all been inserted outside the critical semimajor axis for stability in the 
system. One quarter of the fake transits are not recovered: some fall into data gaps, others into 
noisier parts of the light curve. However, for the case of Kepler- 16, random simulations of the 
Earth-sized planet transits resulted in a 75% recovery rate, using a lower cut criteria compared 
to the one adopted for finding the signals of Kepler- 16b. As seen from Figure HI the inserted 
fakes occupy a clearly-defined cloud, well-separated from the rest, something that be used as a 
merit criterion of what transit signal can be recovered from that particular light curve. This shows 
that the described method works well for finding not only individual transits but, depending on 
circumstances, also for Earth-like transits of quiet stars and even for transits with durations as short 
as only three low-cadence Kepler samples, e.g. the quaternary eclipse of Kepler- 16 represented by 
the red square on Figure HI 



3. Spectra 

3.1. Apache Point Observatory Observations 

KIC 4862625 and Kepler-47 were observed with the 3.5 m telescope at the Apache Point 
Observatory on four occasions between April and July, 2012. Coincidentally, the spectroscopic 
observations of Kepler-47 by Orosz et al. (2012) and ourselves began within 48 hours of each 
other. We used the medium dispersion Dual Imaging Spectrograph (DIS) in its highest resolution 
spectroscopic mode. The pair of B1200/R1200 gratings in combination with the 1.5" slit delivers 
a spectrum with resolution R~3000 and covering simultaneously two 1200A windows centered on 
4500Aand 6500A, respectively. The slit was oriented along the parallactic angle. Each night we 
obtained one or two exposures of each target supplemented by exposures of spectrophotometric and 
radial velocity standards. Nightly we also obtained several HeNeAr lamp spectra for wavelength 
calibration, and we used telluric lines in the observed spectra to correct for offsets due to flexure or 
other instrumental effects. Conditions of all four nights were not strictly photometric. Each target 
was observed for ~900 seconds per night, yielding a peak signal to noise ratio in the continuum of 
20 to 30 per resolution element. 

The data reduction included bias and flat-field correction, aperture extraction, wavelength and 
flux calibration. We compared the APO long-slit spectrum of KIC 4862625 with a library of stellar 
spectra (Pickles 1998). For the R~3000 APO spectra, the FWHM ~45 km s"^ broadening evident 
in the SOPHIE spectra (Section is unresolved. We estimate from the shape of the continuum 
and strengths of particular spectral lines that the best match is spectral type F8 IV. The subgiant 
classification is consistent with the density of star A determined later from the light curve and 
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Fig. 2. — The transit/anti-transit diagram for KIC 4862625. 



- 11 - 



KOI 10020425 dip vs blip 




1.0 1.5 2.0 2.5 3.0 
log[Chisq(const) - Chisq(dlp)] 



3.5 



Fig. 3. — The transit/anti-transit diagram for KIC 10020423, a.k.a. Kepler-47. 
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Fig. 4. — The transit/anti-transit diagram for Kepler- 16 with simulated transits of an Earth-sized 
planet superposed (yellow). 
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Fig. 5. — Planetary transits of KIC 4862625. The normalized and detrended flux centered on 
the five planetary transit events are shown (top to bottom, labeled 1, 5), along with model fits. 
Events 1 to 4 have been offset vertically for clarity. A primary eclipse immediately after egress of 
event 2 has been excised from the plot for clarity. 
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the mass function. The Kepler Input Catalog (KIC) tends to over estimate log(g) for subgiants, 
especially those hotter than about 5400 K, which can lead to underestimates of their radii in the 
KIC, typically by factors of 1.5 to 2 (Brown et al. 2011). For KIC 4862625, the KIC radius 
estimate is 0.806 R©, indeed approximately half the size that we estimate in this work. 



3.2. SOPHIE observations and data reduction 

The two targets were observed at the end of summer 2012 with the SOPHIE spectrograph at 
the 1.93-m telescope of Haute-Provence Observatory, France. The goal was to detect the reflex 
motion of the primary stars due to their secondary components through radial velocity variations. 
SOPHIE (Bouchy et al. 2009) is a fiber-fed, cross-dispersed, environmentally stabilized echelle 
spectrograph dedicated to high-precision radial velocity measurements. For such binary systems 
the amplitudes of variation are expected to be of the order of a few to a few tens km s^^ which is 
well in SOPHIE capabilities despite the faintness of the targets. The data were secured in High- 
Efficiency mode (resolution power i? = 40000) and slow read-out mode of the detector. In order to 
reach signal-to-noise ratio per pixel of the order of 10 at 550 nm, exposure times ranged between 
1200 and 2000 sec for Kepler-47, and between 500 and 900 sec for KIC 4862625. 

The spectra were extracted from the detector images with the SOPHIE pipeline, which in- 
cludes localization of the spectral orders on the 2D-images, optimal order extraction, cosmic- 
ray rejection, wavelength calibration and corrections of flat-field. Then we performed a cross- 
correlation of the extracted spectra with a G2-type numerical mask including more than 3500 
lines, and finally measures the radial velocities from Gaussian fits of the cross-correlation func- 
tions (CCFs) and the associated photon-noise errors, following the method described by Baranne 
et al. (1996) and Pepe et al. (2002). For Kepler-47 and KIC 4862625 respectively, the fuU widths 
at half maximum of those Gaussians are 12 ± 1 kms^^ and 15 ±2 kms^^ and their contrasts 
are 17 ± 4 % and 4 ± 2 % of the continua. One of the observations of KIC 4862625 was made 
at the twilight: the pollution due to the bright sky background was corrected thanks to the refer- 
ence fiber pointed on the sky (e.g. Hebrard et al. 2008). Other exposures were not polluted by 
sky background nor Moon light. In Table |9] the SOPHIE radial velocities are absolute in barycen- 
tric reference, whereas the APO radial velocities are absolute for KIC 4862625 and relative for 
Kepler-47. 

Radial velocity variations in phase with the Kepler ephemeris are clearly detected. Section 
14. 3 1 addresses the orbital parameters derived from the radial velocities in combination with the light 
curve. Section |431 addresses atmospheric parameters of star A derived from the spectra. 
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4. Data Analysis 



4.1. Eclipsing Binary Light Curve 



The light curve of the EB star constrains the relative orbit of the two stars: the orbital period P, 
the center times of primary eclipse Tj and secondary eclipse To, the semi-major axis of the relative 
orbit in units of the radius of star A, a/r^, and the orbital inclination, i. Also from the light curve, 
we derive the relative radii of the two stars, k = rs/rA, and the fraction of the flux in the Kepler 



examine that assumption later in this section. We adopt the period P = 20.000214 d from Prsa et al. 
(2010), although the trend in the eclipse timing variations suggest a period ~2 seconds longer. The 
free parameters of the fit to the light curve are T?, To, a/rA,p, and/g. We compute the fraction of 
light blocked by one star by the other using computer code of Mandel & Agol (2002). Because the 
latter code was designed for planetary transits and models the nearer body as an entirely dark and 
opaque circular disk, we account for the light from the nearer star appropriately by superposition. 

We estimate the limb darkening parameters (mla,"2,a) = (0.243,0.371) and (Mi,g,M2,s) = 
(0.107,0.350) for stars A and B respectively as appropriate for their effective temperatures, grav- 
ities, and metalicities (Table O. For star A, those are estimated from spectroscopy (Section |4^ . 
For star B, we estimate its effective temperature 7^// = 3390 ±50 K, by interpolating a grid of 
model atmosphere spectra (Hauschildt et al. 1995Q) to match the depth of the light curve at sec- 
ondary eclipse; we integrate the spectra over the Kepler bandpass and compare to that of star A, 
accounting also for the relative solid angles, k^. We estimate star B's gravity log(g) = 4.9 in cgs 
units from its mass 0.28 ± 0.034 M© and radius 0.34 ±0.013 Rq. We assume star B's metalicity 
is equal to that of star A. 

We derive the quantity e{cosV)) j Vl — from the phase of the center of the secondary eclipse 
relative to the center of primary eclipse. For the latter constraint, we used the analytic approxima- 
tion for the case of inclination i = 90° (Hilditch 2001, Eq. 5.67), and verified that the difference 
between that approximation and the numerical estimate for i > 87° is negligible (Hilditch Eq. 
5.63). Similarly, the quantity esind) equals the ratio of the difference to the sum of the durations 
of secondary and primary eclipses (Hilditch Eq. 5.69). However, because the eclipse durations 
are much less precisely measured than the centers, we do not explictly constrain esinod, although it 
is weakly constrained implicitly in fitting the light curve. Instead, (0 is measured better using the 
radial velocities (Section l431) r I 




emitted by star B, /g. Nominally, we assume zero "third light," so /a — I — /s; we 



^http://keplergo. arc. nasa.gov/CalibrationResponse.shtml 

^Machine-readable tables are available at |http://svo.cab.in ta-csic.es/theory/db2vo/index.php . 

^We have adopted the equations and viewing geometry of Hilditch's textbook. Figure 2.5. Apparently, prior pub- 
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Given the above constraints from the light curve, fitting the radial velocities depends on only 
three astrophysical free parameters: the systemic velocity y, the velocity semi- amplitude ki, and 
the longitude of periastron co. In our analysis, once co is measured using the radial velocities, the 
eccentricity e and the time of periastron passage tperi are analytically constrained from the light 
curve (Hilditch Eqs. 4.10 and 5.67). Similarly, the uncertainties in the parameters e and tperi flow 
down from the uncertainty in (O determined from the radial velocities. 

The planetary transits of star B, if they occur, are undetectable in the Kepler photometry. 
Because star B contributes only = 0.00124 of the flux of the system, and its relative solid angle 
= 0.0426, the mean surface brightness of star B is 0.029 that of star A, in the Kepler bandpass. 
Given that the planet blocks ~0.1% of the system's light when transiting star A, we predict only 
a 30 ppm planetary transit of star B, i.e. undetectable with Kepler's per-cadence RMS noise of 
222 ppm. Because CB transit durations can be no longer than half of the orbital period of the EB 
(Schneider & Chevreton 1990), or 10 days in this case, and even with an idealized 10-day upper- 
limit on the duration, the 30-ppm transit depth would correspond to ~3a, and typically not even 
that, due to long-term intrinsic variations in the system's total light. 

To estimate eclipse time variations, ETVs (Figure|7]), we empirically determine the best-fitting 
center time of each primary eclipse individually by simply minimizing x^, while adjusting only the 
center time of the model curve. The model is simply the Mandel & Agol (2002) curve that best- 
fits the ensemble of primary eclipses, although simpler, non-physical models such as Gaussians 
or U-shaped curves produce very similar ETVs. We evaluated the statistical significance of the 
measured ETVs by simulation. We superposed model eclipses upon the detrended Kepler data at 
random phases far from the actual primary or secondary eclipses, and produced simulated ETVs 
that are slightly smaller in amplitude than the actual ETVs, implying that the measured ETVs may 
have a contribution from the gas-giant planet but are essentially consistent with noise (Figure |7]). 

Figure [13] and [121 are schematic scale drawings of the EB system. Figure [13] illustrates the 
sizes of the three objects (star A, star B, and the planet) relative to each other and to the barycentric 
orbit of star A, which stretches only a few stellar radii across. The schematics also illustrate the 
positions of the stars at the epochs of each of the five planetary transit events. Because we have 
approximated the transit chord as the diameter of star A for all five planetary transit events, we 
have not illustrated the orbit of the planet in the schematic diagrams. 

In the above analysis, we have assumed that the Kepler photometer records the sum of the light 
from the two stars, star A and star B, and nothing more, i.e. zero "third light." In our spectroscopic 
observing, we took care to inspect images from the acquisition cameras, and noted no stars of any 



lications of CB planets have adopted the opposite viewing geometry, e.g. Fig. 7 of Murray & Correia (2010). This 
affects CO by 180°. 
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significance within the range 1 —2" of KIC 4862625; 2MASS images support this, also. Because 
poor weather thwarted an attempt at adaptive optics imaging of the environs of KIC 4862625, we 
were unable to inspect within the ~1" seeing limit. In general, however, imaging can never prove 
there is zero third light, because any system could be a hierarchical triple star. We investigated 
the effects of assuming a given third-light fraction, /c = 0, 0.1, and 0.2, of the total light. Table 
miists the fc = solution: k = vb/ta = 0.212, a/r^ = 21.95, and i = 87.57°, and that solution is 
used throughout the paper. For /c = 0.1 those parameters are k = rg/rA = 0.228, a/rA = 22.31, 
and / = 87.59°, and for fc = 0.2 they are = re/r^ = 0.251, q/ta = 23.01, and / = 87.66°. Thus, 
compared to the nominal /c = solution, non-zero third light implies larger stellar densities and a 
larger star B relative to star A. 



4.2. Spectral analysis 

The SOPHIE spectra of KIC 4862625 without background pollution were co-added for spec- 
tral analysis. The Ha and H(3 lines were used to determine the effective temperature Tgff = 
6200 ± 150 K. This estimate was made on each line independently in order to check the consistency 
of the results. The spectrum lacks prominent spectral features due to its broad lines combined with 
the low signal-to-noise ratio. This prevents us from carrying out a detailed spectroscopic analysis. 
We could not derive accurate estimate of the surface gravity from the Mg I triplet and the Na I dou- 
blet. The estimation from these lines is log^ ^ 4.0 ± 0.2; it is a typical value for main sequence 
and subgiant stars in that Teff range. We do not find any evidence of Lithium in the spectrum 
nor any sign of chromospheric activity in the Ca ll H and K lines, but the Ha core shows some 
variable emission features. From the width of the CCF we derived vsin/^, = 31 ±2kms^^ and 
[Fe/H]~ -0.15. 

Using these values of Tgff, log^, and [Fe/H], we estimated the mass and radius of the star by 
comparison with a grid of STAREVOL stellar evolution models (A. Palacios, priv. com.; Lagarde 
et al. 2012). We generated a series of Gaussian random realizations of Tgff, [Fe/H] and log^, and 
for each realization we determined the best evolutionary track using a minimization on these 
three parameters. We found = 1 .23 ± O.2OM0, i?* = 1 .70 ± 0.25 R©, and an isochronal age of 
2.613:6 Gyr. 

4.3. Orbital solution of the binaries 

For each targets, the APO and SOPHIE radial velocities were fit simultaneously with a Ke- 
plerian model. The fits are mainly constrained by the SOPHIE data, which are more numerous 
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Fig. 6. — Kepler light curves of KIC 4862625. The normalized and detrended flux of the primary 
and secondary eclipses are shown with respect to the orbital phase of the EB, along with the binary 
star model. The secondary eclipse data have been centered at zero phase for comparison with 
the primary eclipse data. The secondary eclipse data and the residuals (above) have been offset 
vertically for clarity. 
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Fig. 7. — Eclipse timing variationsof KIC 4862625. The observed times of each primary eclipse 
minus the calculated times for a linear ephemeris are shown versus time (an "0-C" diagram). The 
five planetary transit events are indicated by vertical dashed lines. A single primary eclipse is 
missing from the sequence at BJD = 2455567.8. 
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and accurate. The APO radial velocities are much less accurate, but agree with the orbital solution 
derived from the SOPHIE data. 

In addition to the measured radial velocities, the fits take into account the three constraints 
derived from the Kepler photometry: specifically the orbital period P, and the mid-times of tran- 
sit and occultation, Tt and To. The latter parameters strongly constrain e{cos(i)) / \J\—e^ (Section 
14.11) . The radial velocities confirm the orbital eccentricities and allow e and (O to be measured 
individually. In comparison to the radial velocity uncertainties, the uncertainties in the three pa- 
rameters derived from photometry have negligible effects on the final uncertainties of the derived 
orbit parameters. 

The fits to the radial velocites were made using the Levenberg-Marquardt method, and the 
confidence intervals around the best solutions were determined both from variations and Monte 
Carlo simulations. The histograms of the obtained parameters have a single peak and nice Gaussian- 
like appearance. The derived values and uncertainties are reported in Table [2l the best fits are 
over-plotted with the data in Fig. [8l 

The semi-amplitude of the radial velocity variations of KIC 4862625 is ^ = 18.06 ±0.48 
kms^^. With an adopted primary mass = 1.14±O.14M0, this translates in a secondary mass of 
mg = 0.28 ±0.034 M©. Because KIC 4862625 is single-lined spectroscopic binary, star B's mass 
depends on star A's and their uncertainties are similarly coupled. For KIC 4862625 the dispersion 
of the residuals of the fits are 330 ms^^ and 3.4 km s^^ for SOPHIE and APO data, respectively. 
The dispersion is similar or even smaller than the expected uncertainties on the measured radial 
velocities. Uncertainties on the SOPHIE radial velocities might be slightly overestimated. We did 
not reduce them however in order to be conservative. We cannot detect any significant drift in 
addition to the reflex motion due to the binaries. For KIC 4862625 we estimate an upper limit 
±10km s^^ yr^^ for any additional drift. 

Our derived parameters for Kepler-47 agree with those derived by Orosz et al. (2012). The 
semi-amplitude of the radial velocity variations of Kepler-47 is^ = 31.18±0.12 km s^ ^ With an 
adopted primary mass from Orosz et al. (2012), = l.O4±O.O6M0, this translates in a secondary 
mass of niB = 0.357 ±O.O13M0. Because there are so few RV measurements for Kepler-47, their 
dispersion is less than the expected measurement uncertainty. For Kepler-47 the dispersion of the 
residuals of the fits are 25 ms^ and 2.6 km s^ for SOPHIE and APO data, respectively. 



4.4. Doppler Boosting 

The oscillating radial velocity of the primary star of KIC 4862625 is apparent in the Kepler 
photometry. Due to Doppler boosting, when the star is moving toward the Earth, its observed flux 
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Fig. 8.— Radial velocity measurements of KIC 10020423 (Kepler-47, left) and KIC 4862625 
(right) with l-o error bars as a function of time (upper) or orbital phase (lower) together with their 
Keplerian fit and residuals of the fit. The data are from APO (blue, filled diamonds) and SOPHIE 
(red, empty diamonds). 
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increases and when the star is moving away, its observed flux decreases. Figure |9] iUustrates the 
modulation in observed flux as a function of the orbital phase of the EB star. To reduce the effect 
of the rotational modulation on the light curve, we used the mean flux level estimated at each 
point from the sine wave fit at each point (Section [27T1) . We grouped the results in 100 uniformly 
spaced bins in orbital phase. Hence, each point in Figure [9] represents the median of ~ 350 Kepler 
measurements. With an RMS of 222 ppm per original Kepler observation, each median would have 
a formal uncertainty v350 less, or 12 ppm if there were no trends in the light curve. However, 
with the trends, the observed RMS deviation of the medians with respect to the best-fitting boosting 
curve is 19 ppm. 

Due to Doppler boosting, the ratio of observed flux to emitted flux Fq x is 

^ = i-A (1) 

where is the stellar radial velocity, c is the speed of light, and the Doppler boosting factor 
B = 5 + dlnFx/dlnX (Bloemen et al. 2010; Loeb & Gaudi 2003). For a T=6150 K blackbody ap- 
proximation to the star A's spectrum, and a monochromatic approximation to the Kepler bandpass 
ofX = 600 nm, the boosting factor Bbb = 3.99 (Loeb & Gaudi 2003, Eqs. 2 and 3). At a finer level 
of approximation, using a template spectrum for an F8 IV star (Pickles 1998) and the Kepler band- 
pass, we estimate a photon weighted bandpass-integrated boosting factor Bpgiv = 3.73 (Bloemen 
et al. 2010, Eq. 3). Both estimates neglect reddening, but its effect is very small for interstellar 
reddening typical of Kepler stars (Bloemen et al. 2010). Figure [9] compares the Doppler boosting 
effect estimated with B = Bfgjy, to the Kepler photometry of KIC 4862625. The boosting factor 
that best fits the Kepler photometry and the spectroscopic radial velocity curve is 5 = 3.46 ± 0.065, 
i.e. or 93 ± 1.7% of the analytic estimate with B = Bf^jy. The boosted flux from star B is out of 
phase with that of star A, but for simplicity of this analysis we have neglected the tiny contribution 
from star B. 

The capability to measure the radial velocity of an EB star using Kepler data alone could be 
useful and convenient. In principle, CB systems could be "solved" without spectroscopically de- 
termined radial velocities. The phase difference of the stellar eclipses constrains well the quantity 
ecos{(X>) / y/l —e^; so the eccentricity e and the longitude of periastron, (O, are constrained by the 
Kepler photometry. Also, the ratio of the difference in durations (secondary eclipse minus primary 
eclipse) to the sum of the two durations equals esin{(a), to a good approximation that the orbital 
inclination is ~ 90°. With Doppler boosting, the Kepler photometry provides the equivalent of a 
single-lined spectroscopic binary: the radial velocity of the brighter star as a function of orbital 
phase. The latter also constrains e and CO and the phase of periastron passage, along with a measure 
of the radial velocity semiamplitude. For stars of similar brightness, the Doppler boosting effects 
of the two stars will tend to cancel. Of course, the traditional light-curve analysis of the eclipses 
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provides an estimate of the relative brightnesses of the two stars in an EB. For those eclipsing bina- 
ries with one star much brighter than the other, and if that brighter star is photometrically stable or 
at least predictable as in the case of KIC 4862625's rotational modulation, Doppler boosting curves 
from Kepler photometry may provide radial velocities adequate for estimating the mass function 
of the system, and other parameters of the EB. 

In this work, we have demonstrated that the stellar radial velocities can be measured either 
with a spectrometer or a photometer. We give priority to the spectroscopic technique because of 
its well-calibrated heritage. For comparison we measured ki = 16.7 ± 0.5 and (O = 222° ± 2° from 
Doppler boosting in the Kepler light curve. We estimated the maximum-likelihood values of ki 
and (0 from all of the photometry. However, because correlations are apparent in the residuals 
of the Doppler boosting curve, induced by the 2.63-day averaging window during detrending, we 
estimated the uncertainties of ki and (O from subsets of points selected to be independent from 
each other, in steps of phase equal to 2.63 d / 20 d. The photometrically-determined value of (0 
is consistent with that determined from shifts of spectral lines, (O = 220.2° ±3°. Formally, the 
photometrically-determined value for (O is more precisely determined than the spectroscopically- 
determined value. The value of ki is prone to systematic error; it is directly proportional to the 
boost factor B estimated from the overlap integral of the stellar spectrum and the Kepler bandpass. 
The systemic velocity y is indeterminate from Doppler boosting. 



4.5. Stellar Rotation and Star Spots 

From the broadening of the spectral lines and the period of amplitude modulations in the light 
curve, we infer rotation of star A and determine its radius, 

^ P^^vsin/ 

where (|) is a factor of order unity that would account for differential rotation and any systematic 
errors in v sin/*, such as sinz* < 1. With v sin z* = 31 ±2kms^^, p,-ot = 2.6382 ±0.0037 d, and (|) = 
1 the stellar radius Ra = 1 .62 ± 0. IRq, in agreement with the estimates from the spectral analysis 
(Section |43I) and the photodynamical model (Section [6]). With Ra determined from rotation, the 
density of star A implied from the light curve and the mass function from the radial velocity 
semiamplitude, we derive a mass Ma = 1.2±O.2A/0, again in agreement with estimates from the 
other two methods. We also estimated Ma using the predicted velocities of the primary and the 
planet during the five transits from the best- fit simulations. Each of the five events is scaled for 
the respective velocities, such that they have approximately the same width. All five are then fit 
together. The derived parameters of the primary star agree with the above mentioned values within 
their uncertainties. Because M R^, the fractional uncertainty in stellar mass is three times larger 
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Fig. 9. — Doppler boosting of KIC 4862625. The relative, normalized flux change (see text) is 
plotted with respect to the EE's orbital phase. The data binned in 0.01 intervals of phase (filled 
circles) approximately match the flux change estimated from Doppler boosting (red curve) based 
upon the spectrum of the primary star and its spectroscopically-determined radial velocity curve. 
The amplitude of the best-fitting curve (blue) is 93 ± 1.7% that of the estimate (red). 
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than the stellar radius' fractional uncertainty, which is dominated by the uncertainty in vsinz* and 
any bias implicit in (|). 

A curious situation is possible with KIC 4862625. Star A's peak orbital velocity transverse to 
the line of sight is ~18 km s^^ and its rotational velocity is 31 km s^^ (Section 14^21) . Assuming that 
star A's rotation is prograde with respect to the orbital motions of star B and the planet, and the axis 
is in the plane of the sky, i.e. sin i = 1, and that the transit is a central one, the transverse velocity of 
a star spot (or plage) could briefly match that of the planet (~46 km s^^)! Planetary transits could 
exhibit asymmetric shapes due to the nearly matched transverse speeds. KIC 4862625's transit 
events 1, 3, and 5 in particular appear asymmetric, with egress being more abrupt than ingress. 
We are unsure whether the asymmetry is astrophysical or an artifact of the detrending. Because 
they are grazing, the primary eclipses of KIC 4862625 may not exhibit similar asymmetries if 
spots (or plage) are not prevalent near the pole of star A, and in any case the projected rotational 
speeds will be small at the star's poles. Silva (2003, 2008) and Nutzman, Fabrycky, & Fortney 
(2011) have analyzed star-spot induced asymmetries in transit light curves but did not explicitly 
consider the possibility of a star spot over-taking a planet. More typically, e.g. for a 3-day period 
planet transiting a solar-type star rotating every 30 days, the transverse velocity ratio (planet/star) 
is ~100. Gravity darkening is another mechanism for which very rapid stellar rotation can induce 
subtle asymmetries in transit light curves (e.g. Barnes, Linscott, and Shporer 201 1). 



5. Diagnosing a System 

Initial diagnosis of a single transiting planet that orbits a stellar binary can be challenging, 
particularly if orbital period is long and the planet transits only one star. There are many system 
parameters and potentially only a few observational constraints. In some cases the existing data 
will not fully constrain the system. Given this potential complexity, it is useful to understand the 
sequence of analysis steps that build understanding of a newly discovered system. With this in 
mind, we describe the clues we used to diagnose the KIC 4862625 system. 

First, we worried that KIC 4862625 could be an astrophysical false positive. The aperiodicity 
of the five planetary transits disproves the hypothesis of a background eclipsing binary mimicking 
a transiting CB planet. Superficially, the false-positive of a dilute EB mimicking planetary transits 
of a single star (Brown 2003) has an analogy in CB planets, namely that a dilute eclipsing triple 
star in proximity to an EB star could in principle mimic some of the characteristics of a transiting 
CB planet. First, however, the chance proximity on the sky of a double star and a triple star will 
be much rarer than that of a single star and a double star (Lissauer et al. 2012). Second, while one 
could contrive a dilute triple star to mimic the aperiodic centers of a few transit-like events, the 
durations in general would not match also, because for a CB planet, the transit durations depend 
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critically on the characteristics of the EB (Eq. H]). We conclude that a third body orbits the EB star 
KIC 4862625. 

We next analyzed the stellar binary, which has best observational constraints. The Kepler 
light curve folded on the 20 day period of the stellar binary (Figure 16]) has a primary eclipse that 
is 1.3% deep and a secondary eclipse that is 0.1% deep. Primary eclipses are more V-shaped than 
U-shaped, suggesting that a smaller secondary grazes a larger primary. Secondary eclipses have a 
flat bottom with gradual ingress and egress, suggesting that the limb of the primary fully occults 
the secondary. The phase difference between primary and secondary eclipses indicates an orbit 
with significant eccentricity. 

Next we examine high-resolution spectra of the system. Spectral features are roughly similar 
to the Sun, except that the lines are shallow and broad. Cross-correlation with a template yields 
radial velocity shifts consistent with the light curve period of 20 days. Line widths imply vs'mi^ = 
31 km s^^ which is typical for spectral types slightly earlier than the Sun. Fitting the measured 
radial velocities yields eccentricity e = 0.2 and semiamplitude ^ = 18 km s^^ The mass function 
is then 0.008 M©, which means 0.23 < Msec < 0.28 Mgfor 1.0 < Mpn < 1.3 Mq. As expected, 
lines of the low-mass secondary are not detected in the optical spectra. 

Next we examine the five planetary transits that exist in publicly available Kepler data. First 
we checked the data quality flags to verify that all data are valid during the transits. Using the 
PyKE softwarel^ with custom apertures to analyze Kepler target Pixel files, we confirmed that the 
five transit events come from the central target pixels. We measured no significant centroid shifts 
during the transits, which rules out certain astrophysical blending scenarios. We measured transit 
start and end times by fitting U-shaped functions. We then calculated transit midpoint times and 
durations. 

The cadence and duration of the transits suggest that all involve the primary star. The cadence 
of the planetary transits yields an apparent orbital period of approximately 136 days. The time 
differences between successive transits is 136.6, 136.7, 135.9, and 133.2 days. Such regularity 
is not anticipated for a circumbinary system, where a single object transits or eclipses a "moving 
target" (Orosz et al. 2012). The regular cadence of transits is not caused simple commensurability, 
given the 20 and 136 day periods in the system. 

We use the durations of the transits to constrain the parameters of the stellar binary. Transit 
duration depends on chord length and on transverse velocity of the stellar primary relative to the 
circumbinary object (Schneider & Chevreton (1990). During transits, transverse velocity of the 
planet is always fairly similar, but transverse velocity of the occulted primary (dependent of the 



^http://keplergo. arc. nasa.gov/PyKE.shtml 
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mass ratio of the two stars) changes direction and amplitude throughout its orbit. The amplitude 
of the transverse velocity of the primary decreases to a minimum at orbital turning points. For ec- 
centric stellar binaries, amplitude also decreases with increasing binary separation. Finally, transit 
durations are longer when the primary is in the portion of its orbit where transverse velocity of the 
primary and planet are aligned, and shorter when the two velocities have opposite signs. 

With these factors in mind, we now interpret the five observed transits, ordered by increasing 
duration. Figure [121 illustrates the configuration of the stellar binary at the time of each transit. 
Transit 2 has the shortest duration because the transverse velocities are oppositely directed and the 
primary is far from a turning point. Transit 1 has a slightly longer duration because the primary 
is slowing as it approaches a turning point. Transit 3 is near the other turning point, where binary 
separation is larger, further reducing relative velocity. Transit 4 has a relatively long duration 
because the transverse velocities are now aligned, though the primary is near the turning point 
where binary separation is large. Transit 5 has the longest duration because velocities are still 
aligned and the primary is near the turning point where binary separation is small. Even longer 
transit times are possible between points 4 and 5 in Figure[l2l Because transit durations are shortest 
near primary eclipse, the planetary orbit must be prograde relative to the stellar binary. 

To quantify constraints provided by observed transit times and durations, we calculated the 
locations of all bodies in the system as a function of time. Figure [T2l shows the view perpendicular 
to the orbital plane, while Figure [13] shows the view from Earth. The primary has a very similar 
projected position during transits 1 and 5. Because of this near coincidence, the time difference 
between transit 1 and 5 is almost exactly four times the 135.59 day orbital period of the planet. 

Next, we inspect event 4, also a relatively long duration event. At the phase of event 4, the 
primary star is near the opposite turning point of its orbit from events 1 and 5. Note event 2 is 
near primary eclipse, when the two star's relative motions are relatively fast, and the primary star 
is traveling right to left, i.e. in the opposite direction to the planet, if in our scenario the planet is 
traveling in a prograde orbit (left to right in this model). Again, event 2 is deeper than secondary 
eclipse, and event 2's shape is not a square wave, so we consider event 2 is again a transit of the 
planet across the primary star. 

Transverse velocity of the primary is a function of true anomaly, argument of periastron, and 
eccentricity. Using the nomenclature of Hilditch (2001), transverse velocity Vx is given by 

_ M2 ^ I iTlGMbin ^ 1 /3 (e sin (D + sin(9 + CO) ) 

^mJ^ Pun ^ (1 -.2)1/2 ' 

where A/^;>j = Mpri -|-A/sec, 9 is true anomaly, (O is the argument of pericenter and e is the eccen- 
tricity of the binary star. 



Using the formalism of Schneider & Chevreton (1990), and assuming a circular orbit for the 
planet, we obtain the following expression for transit duration tdm 
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X = (e sin (0 + sin (0 + 0))) 



In Eq. m we use the mass function /(m) obtained from radial velocities to substitute secondary 
mass (M2) for binary mass (Mbin). The mean period of the planet (Pp) is known from the cadence of 
transit times. Assuming transit chord (Re) is the same for every transit, we use amoeba/curvefit to 
fit for the coefficients A and B (C is known to the precision of Pp). The inferred value of parameter 
B in Eq. |4] constrains stellar radius. Long duration transits will deviate slightly from Eq. |4l if 
velocity of the stellar primary changes significantly during the transit. We caution that very small 
changes in observed transit durations can have significant effects on parameters derived from A 
and 5. 

To validate Eq. |4l we fitted observed transit durations for Kepler-47b, using transit parameters 
in Orosz et al. (2012). The fit shown in Figure [TOlyields Mhin = 1-35 Moand Rp^ = 0.87 Rq, which 
agrees fairly well with 1.4 Mq and 0.96 Rq deduced by Orosz et al. (2012). Figure [TT] shows an 
analogous fit for KIC 4862625. We obtain an Mfe,„ = 1.34 Mq and Rpn = 1.8 Rq. The analytic 
curve indicates that future transits may be as long as a day! To assess uncertainties, we created 
simulated observations for a mock system like KIC 4862625 (Model 1 in Section |6l). We integrated 
the mock system for 9 planetary orbits, using a time-step of 14.4 min. Blue diamonds in Figure [TT] 
show the result. As with Kepler-47b, the inferred binary mass of 1.35 Mqj is ~ 8% larger than the 
expected value of 1.26 Mq. The chord length calculated from A is 1.9 Rq, compared to the input 
value of 1.8 Rq. The fits to the planetary transit durations using the analytic functions provide a 
good starting point for the more refined fits presented in Section |6l 

In the existing data, depths of the five observed transits appear to change with time. Perhaps 
this is simply an artifact of our detrending procedure. Alternatively, star spots that perturb the 
light curve may also affect apparent transit depth. Finally, changes in transit depth could be a 
manifestation of a planetary orbit with an inclination smaller than 90 degrees. 
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Numerically integrating an inclined planetary orbit forward in time, the transits fade away, 
cease for a period, and then return after a large number of orbits. The larger the inclination, the 
faster the evolution. Given that we have only observed 5 transits so far, we assume that the chord 
lengths, transit depths, and planetary inclination are constant. 

For completeness, we note that the observed dependence of the durations of the five transits 
on the phase of the binary rule out a retrograde orbit for the planet. Such an orbit should exhibit an 
opposite trend in a duration versus phase diagram, namely short durations near secondary eclipse 
and longer ones close to primary eclipses. 



6. Planetary Transits 
6.1. Methodology 

Combined light curve and radial velocity measurements of a circumbinary planet can be char- 
acterized by 17 parameters: five orbital elements of the binary (abm, ^hin, w, Q., i), six osculating 
orbital elements of the planet (as, 63, 23, (O3, ^^3, (|)o), three masses (Mpn, Msec, ms), and three 
radii (^pri, ^sec, fi)- Exhaustively searching a space with 17 dimensions and nonlinear parameters 
is impractical, so we make simplifying assumptions that allow us to obtain a reasonable solution. 
As discussed in Section RV, radial velocity measurements yield e, (O, and the stellar mass ratio. 
The precise Kepler light curve yields the ratio of stellar radii, which depends on derived impact 
parameters for eclipses that may be grazing. We assume the five observed transits all occult the 
primary star, so radius of the secondary star drops out of the system. For simplicity, we assume 
that the planet has negligible mass and orbits in the plane of the sky (ii, = 90 deg). The former 
assumption is based on the measured ETVs, discussed in Section ETV and expanded below. The 
latter relies on the hypothesis that the inclination of the planet should be certainly larger than that 
of the binary for otherwise it would not be seen in transit. 

With these approximations, the four remaining free parameters are M^^, a^,, e^, and (O3. Find- 
ing optimal values for these four parameters still requires a very fine numerical grid. For example, 
a 1% change in mass of the primary star, or in the semi-major axis of the planet can dramatically 
change the dynamical evolution of the planet, affecting arrival times of observed transits by hours 
or even days. Arrival times would also be affected if the planet is massive, rather than a massless 
test particle. 

Transit durations (shown in Table |4] have formal uncertainties that do not necessarily account 
for astrophysical variations in the light curve. To assess how sensitive binary mass is to individual 
transit durations, we generated and analyzed a set of perturbed observations. We scrambled the five 
observed transit durations and added a normally distributed perturbation with a standard deviation 
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Fig. 11. — Similar to the previous figure but for KIC 4862625. The blue diamonds represent the 
predicted transits from the dynamical model (Model 1, Section [6]). The fit is over the data points, 
not over the simulated durations, which are shorter than the observed ones by up to an hour (see 
Table HI). The uncertainties in the measured durations are smaller than the size of the data points. 
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Fig. 12. — Scaled view of the binary system of KIC 4862625 viewed perpendicular the orbital 
plane. The inner ellipse shows motion of the primary, while the outer ellipse shows the secondary. 
Numbers along the orbits indicate times of the three transits. The diagram is to scale. 
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Fig. 13. — Schematic scale diagram of KIC 4862625. The diagram is to scale, with solar radii 
indicated on the axes. The positions of star A in its orbit (black ellipse) are plotted at the times of 
the five planetary transits (labeled, black circles centered on black + symbols). The corresponding 
positions 0.25-day before (blue) and after (red) are also shown, to illustrate the projected velocity of 
the star. Star B is shown to scale (dashed circle, arbitrarily positioned below star A at event 2). The 
planet also is shown to scale (smallest circles); again, blue and red correspond to positions — 0.25-d 
and -|-0.25-d with respect to the position shown in black, which has been arbitrarily centered at the 
origin. 
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of 1 hour. We then used Eq. |4]to evaluate the binary mass, obtaining the results shown in Figure [141 
Except for a few realizations with unrealistically large perturbations, binary mass is concentrated, 
with 68% of the values in the range 0.8 to 1.7 Mq. In subsequent analysis, we constrain binary 
mass to be in this range. 

A simple Keplerian solution cannot reproduce observed planetary transit times because the 
central potential varies as stars in the binary orbit each other. Instead we use a three-body numerical 
integrator with hierarchical Jacobian coordinates, as will be further described in Section |7] We 
specify initial conditions (and reference epoch) with respect to the midpoint of the first transit, an 
arbitrary yet convenient definition. The coordinates of the three bodies are defined with respect 
to the barycenter of system, which is the barycenter of the stellar binary when planet mass is 
negligible. We use time-dependent, osculating Keplerian elements to calculate a^, £3, 13, and (03 
(Doyle et al. (2011)). We use the numerical integrator to compute spatial coordinates for the two 
stars and the planet. Transits intervals occur when the projected distance between the primary star 
and planet centers is less than the sum of their radii. Note that parameters cited in this section are 
instantaneous values for the reference epoch. 

First, we used a coarse time step (one third the shortest transit time) to calculate dynamical 
solutions and corresponding transit midpoint times for a grid of parameter values. The grid step 
sizes were 0.01 Mq for Mpri, 1 Rq for 03, and 5 deg for (1)3. Changes of 1 Rq in aj, cause significant 
changes in transit midpoint times. Parameter ranges in the coarse grid were 0.8 to 1 .7 Mq for Mp^, 
to 0.2 for 63, and to 360 deg for (03. Next, we used a fine time step (0.01 days) to calculate 
dynamical solutions for a much finer grid of parameters around the best solutions in the coarse 
grid. The fine grid includes ri as a fifth parameter and calculates transit durations in addition to 
transit midpoint times. Comparing computed and observed midpoint times and durations for the 
five transits yields a set of plausible models. 



6.2. Results 

Even searching in a restricted parameter space, we found multiple solutions that are consistent 
with the five planetary transit observations. We adopted a goodness of fit metric that is the root 
mean square (RMS) of observed minus calculated (O— C) midpoint times for transits 2 through 5. 
The midpoint time for transit 1 always has zero residual because it defines the zero point of time 
in our calculations. Figure \T5\ shows goodness of fit for every dynamical model in the grid that 
matches observed transit midpoint times to better than a RMS of ~ 7 hours. Each point in the 
figure represents a unique combination of Mp^, aj, 63, and (03. Despite a uniformly sampled grid, 
the number of solutions better than RMS of ~ 7 hours decreases with increasing Mpri, suggesting 
that lower Mp^ is more likely. We note, however, that even though almost half of all solutions fall 
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within Mpribetween 0.8 M0and 0.9 Mq, most of them have high RMS and are not very likely. 

Within the limitations of our grid space search, three solutions stand out from the rest as 
having an RMS O— C midpoint times of events 2 through 5 to within 1 hour. Granted we may have 
fallen prey to a local minimum, the three models, labeled Ml, M2, and M3 in Figure[T5l reproduce 
the observed planetary transits equivalently well, with a RMS of ~ 1 hour. Further supporting 
the existence of a solution around this particular parameter island stems from them having very 
similar parameters - all three indicate a primary star of Mpyi ~ IM© (within lo of the estimates 
from the mass function and the spectra) and a planetary orbit with a semi-major axis a3 ~ 0.57 
AU, and ~ 0.1. Combined with observational constraints described in previous sections, this 
implies Msec ~ O.26M0, Rp^ ~ 1.6Rq (consistent with our estimates from the spectra and the 
stellar rotation), 7?sec ~ 0.34^?©, and ^3 ~ 138 days. Table[3]gives the best-fit parameters for each 
solution. The models reproduce well the observed phase-dependance of the durations, with event 
2 predicted to be shortest and event 5 - longest. The five observed transits are not sufficient to rule 
out any of these plausible solutions. Other equivalent or even better solutions may exist, given the 
coarseness of our initial grid search. We were unable to match observations with a circular orbit 
(^3 = 0), which is not surprising given that orbital elements evolve continuously. An orbit that is 
initially circular with change with time, especially when the planet is relatively close to the stellar 
binary (e.g., KIC 4862625). 

To estimate the uncertainty of the derived Mpi-i, we examined the distribution of solutions 
as a function of Mpj-i that have an RMS of less than 3 hours (~ 2a above the best-fit value). Of 
the 499 solutions on Figure [H only 167 comply with this merit criteria, of which 63 fall within 
Mpri between 0.9Mq and I.IMq. The distribution is double-peaked with a taller one centered on 
Mpri ~ O.8M0 and another, smaller peak near Mpri ~ IMq. As mentioned above, the low-mass 
solutions have significant RMS and, as discussed in Section 1421 we expect Mp^ to be ~ IMq. 
Thus, we focus on the Mp^ ~ IM© peak, where half of the solutions fall within a range of Mp^ 
between O.9M0 and 1 . IMq, the uncertainty value we adopt under the low-number statistic situation 
we face. The same approach help us constrain the semi-major axis of the planet and the eccentricity 
of the orbit but not the argument of periastron where the best-fit values were uniformly spread. We 
report the values in Table [3l 

After global minimization, individual transits should be inspected for significant inconsisten- 
cies between observed and model light curves. To illustrate this. Figures [T6l [TTl and [T8] compare 
our three best models with the five observed transits. Dashed lines indicate observed times of first 
contact, transit midpoint, and fourth contact. Agreement is very good for all three models, but 
there are subtle differences - the solutions oscillate around the observed transits. For example, for 
Model 1 (blue) transit 2 is late, transit 3 is early, and transit 4 is on time. In contrast, for Model 2 
(red) transit 2 is on time, transit 3 is on time, and transit 4 is late. It is not obvious which model 
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is better. The fit is limited by uncertainties in the observed transit midpoint time and resolution of 
the dynamical models. 

To study whether future transits can help break this degeneracy, we continued integrating the 
three models into the future. Table |4] gives the results, listing midpoint times and durations for 
future transits 6 through 9. By transit 9 on day 6157 (BJD-2450000), Models 1 and 2 predict 
transits times that differ by ~ 3.5 hours, which is certainly detectable. Coincidentally, transit 9 
will have a long duration and be near the peak of the analytic function described in Section [5j 
yielding an improved stellar mass constraint. Events 6, 7, and 8 will have short duration and be 
in the "valley" of the duration curve (Figure [TT1) . providing limited additional information about 
stellar masses. We note that Figure [15] may be rearranged with the addition of future transits such 
that our best-fit solutions fall out of grace only to be replaced by other - for example, nearly 50 
more solutions are currently within ~ 2o of the best RMS. 

We also investigated how non-zero planet mass (1 to 10 Mj) affects transit midpoint time and 
the precise cadence of stellar eclipses by perturbing Model 1. Using the formalism of Borkovits et 
al. (2012), we estimate that light travel time effects due to a 10 Mj planet would produce timing 
changes of only a few seconds. On the other hand, dominant three-body dynamical interactions 
would produce an amplitude of ~ 5 minutes, which is twice the observed value (Figure |7]). Thus, 
the mass of a coplanar planet should not exceed ~ 5 Mj. We rule out a more massive body on a 
highly inclined orbit because we observe planetary transits. Figure [19] shows how predicted transit 
midpoint times change as planet mass increases from to 5 Mj. A large planet tends to arrive 
earlier for events 2 and 3 and late for 4 and 5, where the deviations are more pronounced for the 
latter two. We tested several other models and found they all show similar behavior, indicating that 
while Saturn-mass planets and smaller can be safely approximated as test particles for dynamical 
purposes, exo-Jupiters' masses need to be properly accounted for in numerical integrations. 



6.3. Search for additional transits 

We visually examined the light curve of KIC 4862625 for evidence of additional transits. 
Despite significant variability, there are a few features that are reminiscent of shallow planetary 
transits. For these features, data flags are normal and centroids have no significant shifts. Some 
of these extra features occur close in time to the five main transits, most notably around day 5615 
(BJD-2450000) when the dynamical model predicts the planet is near the secondary star, but also 
near days 5045 and 5069, where the two are distinctly apart. The small size and faintness of the 
secondary, however, means we do not expect to see a quaternary transit (see Section [43]) . Others 
(near days 5258, 5269, 5298, 5312, 5415; BJD-2450000) cannot be easily associated with the 
circumbinary planet KIC 4862625b. 
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Fig. 14. — Mass of the binary (M^ derived using the results from the RV measurements and 

applying Eq. |4]to a thousand iterations of the five planetary transit durations, scrambled and with 
added normalized noise with o = 1 hour. 68% of the points fall within binary mass between 0.8 
and 1.7 Mq, providing a manageable search space for the primary mass in our dynamical model of 
KIC 4862625. 
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Fig. 15. — Best-fit dynamical models of KIC 4862625 predicting the midpoint times of transits 
2 through 5 within an RMS of ~ 7 hours. Models with different parameters produce similar 
solutions, all predicting planetary transits that deviate by only a few data points from the observed 
values, a small error margin on the orbital scale of the planet. The decrease in the number of "hits" 
as the mass of the primary increases from 0.8 to 1.6 is not a systematic effect but is in fact real - 
there are many more good "hits" for smaller Mprim- 
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Fig. 16. — The three best-fit dynamical models of KIC 4862625 predicting the midpoint times 
for transits 1 and 2 to within two hours. The dashed vertical lines define the time of first contact, 
midtransit time and time of fourth contact, from left to right respectively. All three models are 
consistent with a primary star of 1.6i?0 radius. While the three models predict similar times of 
transits for events 2 through 5, the differences between them are subtle - the solutions "oscillate" 
from one event to another. Future events will allow us to further constrain the characteristics of the 
system. The parameters of the system for the three models are outlined in Tabled 
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Fig. 17.— 



The same as the previous figure but for the third and fourth transits of KIC 4862625. 
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Fig. 18. — The same as the previous figure but for the fifth transits of KIC 4862625. 
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Fig. 19. — The effect of a non-zero planetary mass on the predicted times of transit of KIC 
4862625. For the same parameters of the system (Model 1), a increasingly more massive planet 
introduces significant deviations relative to the best solutions for a massless planet. While it is 
unlikely that the mass of the planet is larger than a few Mj, it adds a significant complication to the 
best-fit models, again indicating the non-unique nature of the solutions. 
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7. Dynamical analysis and orbit stability 



We have carried out a dynamical analysis of the KIC4862625 system within the framework 
of the three-body problem. Such a dynamical system is well known to possibly exhibit complex 
dynamical behavior. In particular we have carried out a dynamical analysis of the planet around 
the binary pair in order to detect chaotic regions, often associated with mean-motion resonances 
(MMRs), in the orbital parameter space of the planet. 

We have applied the MEGNo|zl factor dcincotta & SimdboOlallbl i lcincotta. Giordano & Simo 



20031). This numerical technique is efficient in distinguishing between chaotic and quasi-periodic 
and has found widespread appli cation within dyna r tiical astronomy an d the dynamics of multi-body 
extrasolar planets dGozdziewski et al.ll200l[l2008l : lHinse et alboiol) . 



Chaotic orbits are usually (but not always) attributed to unstable orbits. For a quasi-periodic 
time evolution of the system the dynamics is regular and characterized by only a few funda- 
mental frequencies often associated with stable orbits. However, in order to be precise, quasi- 
periodic/stable orbits can only be quoted as stable up to the considered integration time. Knowl- 
edge on the subsequent dynamical evolution of the system is still hidden to the experimenter. In 
this work we have experimented with various integration length and found an integration time scale 
that is long enough to detect the most important mean-motion resonances close to the osculating 
orbital elements of the transiting planet KIC4862625b. 

For reasons of optimization and execution spee d we applied the rt i ost recent implementatio n 
of the MEGNO technique using Poincare variables (IGozdziewskill2003l : lGozdziewski et al.ll2008h . 
In this reference frame the coordinates of each body are expressed in an astrocentric frame with the 
most massive object at rest and velocities are relative to the systems barycenter. We use geometric 
Jacobian coordinates as our initial conditions where the planets orbital elements are relative to 
the barycenter of the binary system. The Jacobian cartesian coordinates are then transformed to 
Poincare elements. T he numerical integration o f the equations of motion and the corresponding 
variational equatio ns ( Mikkola & Innanen 1999 ) are based on the ODEx|§ integration algorithm 
jHaireretal.lll993l) . 

To compute MEGNQ maps we used the MECHANIcl^ software ( Slonina. Gozdziewski & Migaszewski 



2012l : ISlonina et al.ll2012l) which utilizes a Message Passing Interface (MPI) based framework that 
implements massive single serial dynamical problems for parallel execution in a multi-processor 



^Mean Exponential Growth of Nearby Orbits 
^http://www.unige.ch/math/folks/hairer/ 
^http://git. astri.umk.pl/projects/mechanic. 
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initial semi-major axis (AU) of planet 

Fig. 20. — Upper panel: MEGNO map of KIC 4862625b. The mass of primary and secondary 
component were set to 1.0, 0.26 Mq, respectively. The binary semi-major axis and eccentricity 
was set to 0.156 AU and 0.2. The planet was treated as a test particle. The cross-hair marks the 
best-fit osculating elements (a,e) = (0.57 AU, 0.1) of the planet. Arrows indicate the location of 
mean-motion resonances. Yellow (or light) color denotes chaotic (possibly unstable) orbits and 
blue (dark) color denotes quasi-periodic orbits with | (7) — 2.0| ~ 0.001 . Lower panel: Same as the 
upper panel but now zooming into a narrower region of (a,e)-space of the planet. See electronic 
version for colors. 
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computing environment. Usually we allocated 60 CPUs to compute a single dynamical MEGNO 
map with typical resolution of 500 x 350 initial conditions. Typical integration time of a single 
orbit covered about 5x10^ binary orbits (~ 2700 years). 

When investigating a dynamical system the MEGNO factor (usually denoted as {¥)) has the 
following mathematical properties. In general, MEGNO has the parameterization {¥) (t) = at + ^ 
(see references above). For a quasi-periodic chaotic initial condition, we have a ~ and P ~ 2.0 
for ? — )■ oo asymptotically. If the orbit is chaotic, then (7 ) — )• X? /2 for ? — )• oo. Here X is the maximum 
Lyapunov exponent of the considered initial condition. In practice, we terminate a given numerical 
integration of a chaotic orbit when (Y) > 12. In our MEGNO maps we color code chaotic initial 
conditions with yellow and quasi-periodic with blue color. 

In Fig.[20]we show a MEGNO map displaying the dynamics of a planet in semi-major axis (a) 
and eccentricity (e) space. The planet was assumed (justified by the planet displaying transit-like 
signals) to be co-planar with the binary orbit and all planetary orbital angles were set initially to 
zero. We remind the reader that the shown elements are referred to a Jacobian reference system 
with the planets semi-major axis being measured relative to the binary barycenter. The secondary 
binary component was also started with all orbital angles initially set to zero. The mass of the 
two binary components were set to 1.0 and 0.26 M© for the primary and secondary component 
respectively. Referring to Fig. [20] we identify the location of several mean-motion resonances in 
the system appearing as vertical chaotic columns at constant semi-major axis. The derived orbital 
elements of the circumbinary planet have (a, e) — (0.57, 0. 1) placing the planet in a quasi-periodic 
region between the 6: 1 and 7: 1 mean-motion resonance. We show a zoom-in plot in the lower panel 
of Fig. [20l For planetary eccentricities larger than e ^ 0.25 the planets orbital pericenter distance 
{q = a{l — e)) starts to become comparable to the binary separation resulting in close encounters 
and hence strongly chaotic and most possibly unstable behavior. We detected collisions of the test 
mass with one of the binary components and/or ejection of the planet from the system by carrying 
out single orbit integrations of the planetary orbit for eccentricities around 0.25 and higher. 

We investigated the significance of the mass parameter of all involved objects on the overall 
dynamics. We recomputed the map shown in Fig. [20] for various masses of the planet by consid- 
ering 1 Mgarih-, 10 Mgarth-, 1 Mj^p, 10 M jup and 50 Mjiip masses and found no significant qualitative 
difference in the global dynamics of Fig. [20| From comparison we found that all maps resemble 
each other more or less qualitatively. In a similar mass parameter survey study we also considered 
different mass combinations of the binary components that seemed most plausible (see previous 
sections) and conclude that only small qualitative changes were observed for small changes in the 
components mass. However, the most significant changes were detected to be close to the location 
of mean-motion resonances. We note that the timescales of precession frequencies of some of the 
orbital elements might change significantly for different mass parameters. In addition, we also 
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note that the ETVs can also change significantly for different masses as discussed earlier. 

We conclude this section by noting that the best-fit orbital parameters locates the transiting 
planet KIC4862625b in a quasi-periodic region in (a, e) -space rendering our solution to be plausi- 
ble from a dynamical point of view. We point out that our work suggest that the planet is relatively 
close to the 7: 1 mean-motion resonance with the binary. We speculate that future observations will 
reveal whether the planet is or is not part of a resonant configuration. 



8. Discussion 

A single transit event indicates the presence of a third body in an eclipsing binary system. 
Two transits can define an approximate orbital period for the third body, albeit with uncertainties 
associated with aliasing and data gaps, or orbital inclination that can prevent a planet from transit- 
ing the moving targets of the two stars on each passage of the planet. Even with a few planetary 
transits observed, unraveling the system configuration can be challenging. 

In general order of increasing challenge, the system configurations for circumbinary planets 
can be solved for double-lined, double-dipped circumbinaries (planet transits are detected across 
both stars, e.g. Kepler 34 and 35), followed by double-lined, single-dipped circumbinaries (plane- 
tary transits are detected across only one of the stars) and the single-lined circumbinaries, double- 
dipped and single-dipped. For double-lined eclipsing binary stars, spectra provide the individual 
stars' masses; for single-lined eclipsing binaries, spectra provide only the mass function that re- 
lates the two masses to each other, not the two masses individually. Transits of a third body can 
break the latter degeneracy inherent to single-lined systems. Single-lined, double-dipped systems 
wherein transits across both stars occur during a single binary orbit, allow excellent constraints 
on the masses of the binary stars, as anticipated by Schneider & Chevreton (1990) and demon- 
strated for Kepler-16 by Doyle et al. (201 1). Whether a binary is single- or double-lined depends 
on observational capabilities; for example, Kepler-16 was originally a single-lined double-dipped 
system; it has since changed its category to a double-lined system due to very high sensitivity 
spectra (Bender et al. 2012). Likewise, whether a system is single- or double-dipped also depends 
on the observational capabilities: even if a planet transits both stars, we may not be able to discern 
the transits of the fainter star, e.g. KIC 4862625. More challenging than the systems that are either 
double-lined or double-dipped are the single-lined and single-dipped circumbinary systems such 
as Kepler-38, Kepler-47 and KIC 4862625. A large number of planetary transits, as in the case of 
Kepler-47b, helps photodynamical modeling to constrain system parameters via Eq. |4](cf. Figure 
[TOl) . However, if only a few transits occur, e.g. Kepler-38 and KIC 4862625, there may be doubt 
as to the uniqueness of a solution with a large number of system parameters. 
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From the broadening of the spectral lines and the period of amplitude modulations in the light 
curve, we infer rotation of star A and determine its radius, which in combination with log{g) from 
stellar spectra or the density of star A from the stellar eclipses, either one indicates a Solar-mass 
primary star. Small- or undetectable deviations from a linear ephemeris in the primary eclipse 
times prove the planet is of planetary, not stellar, mass. 

The five planetary transits further constrain the parameters of the two stars, as their center 
times and durations depend on the ratio of the stellar masses and on the transiting chord length 
(Eq. 131). The latter breaks the degeneracy inherent to the mass function of the single-lined spec- 
tra. Thereby we constrain the individual stellar masses of of KIC 4862625 and Star A's radius, 
assuming central planetary transits. Similarly, we confirm the parameters of the Kepler-47 system 
reported by Orosz et al. (2012). 

To solve the dynamical nature of the KIC 4862625 system we did a minimization search over 
a grid space defined by mi, aj, es, and (03, using a three-body numerical integrator in Jacobian 
coordinates. Taking the first event as a reference point, we found a set of best-fit solutions, defined 
by the system parameters of the four-dimensional grid space listed above, that predict the midtran- 
sit times and durations of the subsequent four events to within an hour. The simulations match the 
observed events well but we caution that the combination of fixed grid resolution, triaged param- 
eters space, and the small number of transits limit our ability to choose one of the best fit models 
over another. Observations of a few additional transits will differentiate our models, because the 
optimal solutions diverge in their predictions for the center times of the planetary transits. 

To detect chaotic solutions in the parameter space, we studied the long term stability of the 
system using an extensive set of numerical simulations, applying the MEGNO factor. We tested 
systems with different planetary masses, between 1 Mearth and 50 Mjup, to evaluate the changes 
in the dynamical behavior of the three bodies. We do not detect significant difference outside 
mean motion resonances; a planet starting near a mean-motion resonance, however, exhibits erratic 
behavior. The ratio of binary star period to the period of the giant planet is, however, not an integer 
value, giving us confidence that the planet is not on a chaotic orbit. Its orbit is near but beyond 
the instability region; the ratio of planetary to binary semi-major axis is ~ 3.6, compared to ~ 2.8 
of Holman & Wiegert (1999), which is similar to the other Kepler planets and in agreement with 
theoretical predictions that such configurations should be common. 



9. Conclusions 

We report the discovery and characterization of a gas giant r = 0.49 ± 0.01 8 rjupiter circumbi- 
nary planet transiting the KIC 4862625 eclipsing binary system. The planet revolves around the 
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two stars every ~ 138 days and transits the 1 . 14 ± 0. 14 Mq and 1 .59 ± 0.06 Rq F8 IV subgiant pri- 
mary. Dynamical models indicate a slightly smaller mass for the primary star (IMq) but a similar 
size (1.6i?0), and a 0.56 AU, e = 0.1 orbit for the planet. We describe a simple, semi-automatic 
procedure specifically designed to search for aperiodic transits in the light curves of binary stars. 
After identifying the transiting signals, we obtained spectra of the two pairs of binary stars, from 
which we determined the mass function, eccentricity and argument of periapsis for each pair. 

We also describe the independent discovery of the Kepler-47bc circumbinary planets by the 
same method. We discontinued our characterization of that system soon after it was reported by 
Orosz et al. (2012). Our truncated characterization confirms the parameters they reported. 

We coin a phrase to describe circumbinary planetary systems: if planetary transits are detected 
for only one star, the system is "single-dipped," and for both stars, "double-dipped." We discuss 
the relative merits of single- or double-lined and single- or double-dipped circumbinary systems. 

Periodic variations in the radial velocity of star A measured by Doppler boosting compare 
favorably with those obtained by the traditional spectroscopic methods (Figure|9l). The example of 
KIC 4862625 demonstrates the potential of the Doppler-boosting technique for reconnaissance of 
eclipsing binary stars prior to, or in lieu of, obtaining high-resolution spectra of them. 

The family portrait of circumbinary planets discovered by the Kepler mission is filling up 
quickly, with now seven planets reported in less than a year since the first was reported by Doyle 
et al. (2011). KIC 4862625 joins its peers Kepler 16b, 34b and 35b in the category of gas giant 
planets and, like Kepler 38b, orbits a binary system that includes an evolved primary star. With the 
continued operation of the Kepler mission and its exquisite-quality data, we expect the discovery 
of circumbinary planets to continue and give us a deeper insight into these exciting systems. 
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Table 1: Measured radial velocities. 
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KIC 4862625 
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4.1 
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4.1 
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56097.8526 
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56111.6009 


-41.0 
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-9.13 


0.45 


OHP193/SOPHIE 


56162.5848 


-9.42 


0.64 


OHP193/SOPHIE 


56175.3481 


-28.37 


0.76 


OHP193/SOPHIE 


56176.3531 


-21.09 


0.62 


OHP193/SOPHIE 


56 177.2970^ 


-15.41 


0.50 


OHP193/SOPHIE 


56185.5063 


-17.26 


0.49 


OHP193/SOPHIE 


56188.4707 


-29.40 


0.64 
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56193.5113 


-43.07 


1.03 


OHP193/SOPHIE 


KIC 10020423 (Kepler-47) 
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-65.9 


4.0 
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56093.8661 


0.0 


4.5 


APO 


56097.8526 


-59.7 


5.1 


APO 


56159.5866 


26.67 


0.08 
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56 160.5787 


34.88 


0.13 
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24.22 


0.08 
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t: measurement corrected for sky background pollution. 
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Table 2: Parameters of the Binary Star Systems. 
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3 
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AU 
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i 
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Radius of Star A 
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Table 3: Best- fit Parameters for Three Dynamical Models. 



Binary Stars 


Model 1 


Model 2 


Model 3 


Uncertainty 


Mass of Primary Star [Mq] 


1.0 


1.02 


1.03 


0.1 


Mass of Secondary Star [Mq] 


0.2614 


0.2646 


0.2662 


0.03 


Radius of Primary Star [Rq] 


1.6 


1.6 


1.6 




Radius of Secondary Star [/?©] 


0.35 


0.35 


0.35 




Semimajor Axis [AU] 


0.1558 
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0.1572 


0.003 


Circumbinary Planet 


Orbital Period [days] 
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Semimajor Axis [AU] 
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90 


90 





Event # Center a Depth' a Duration a Center Duration Center Duration Center Duration O 

(BJD-245()()()()) (Center) [ppm| (Deplli) |days| (Duration) (Modellj (Modellj (Model2) (Model2j (ModelSj (Model3j 

- O 

Observed ^ 
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2 5207.4077 0.011 631 90 0.5125 0.021 5207.43 0.51 5207.42 0.50 5207.42 0.50 ^ 

3 5344.1308 0.012 914 83 0.6184 0.023 5344.10 0.55 5344.12 0.55 5344.11 0.55 1^ 

4 5480.004 0.015 1042 78 0.7513 0.029 5480.01 0.76 5480.04 0.74 5480.02 0.74 g_ 

5 5613.2329 0.013 1043 67 0.891 0.032 5613.24 0.84 5613.21 0.82 5613.19 0.82 k 

Predicted q. 

6 — — — — — — 5749.08 (X51 5748.98 05 5748.97 05 ^ 

7 — _____ 5885.71 051 5885.57 0.5 5885.55 05 g_ 

8 — _____ 6022.20 0.6 6022.05 0.6 6022.02 0.6 ^ 

9 — — — — — — 6157.40 1.0 6157.24 1.0 6157.19 1.0 



